Astronomy & Astrophysics manuscript no. Remass 


© ESO 2010 


August 19, 2010 





Mass reconstruction by gravitational shear and flexion 

Xinzhong Er 1,2 , Guoliang Li 1 ' 3 , and Peter Schneider 1 



o 

(N 

< 
OO 



o 

U 

6 

CO 



> 
oo 
oo 
o 
cn 

oo 

o 
o 



1 Argelander-Institut fur Astronomie, Universitat Bonn, Auf dem Htigel 71, D-53121 Bonn, Germany e-mail: 
xerSastro . uni-bonn . de 

2 International Max Planck Research School (IMPRS) for Astronomy and Astrophysics Auf dem Hugel 69, D-53121 
Bonn, Germany 

3 Department of Physics, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907, USA 



Received -; accepted- 



ABSTRACT 



Aims. Galaxy clusters are considered as excellent probes for cosmology. For that purpose, their mass needs to be 
measured and their structural properties needs to be understood. 

Methods. We propose a method for galaxy cluster mass reconstruction which combines information from strong lensing, 
weak lensing shear and flexion. We extend the weak lensing analysis to the inner parts of the cluster and, in particular, 
improve the resolution of substructure. 

Results. We use simulations to show that the method recovers the mass density profiles of the cluster. We find that 
the weak lensing flexion is sensitive to substructure. After combining the flexion data into the joint weak and strong 
lensing analysis, we can resolve the cluster properties with substructures. 
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1. Introduction 

Galaxy clusters are the most massive virialized structures 
in the universe, and are considered as valuable probes for 
cosmological parameters. From numerical simulations one 
can identify collapsed structures (called halos) and study 
their properties. In ACDM simulation, the density profile 
of halos has a universal form, closely following the NFW 
profile (Navarro et al. 1996). Moreover, numerical simula- 
tions predict a rich population of sub-clumps within the 
virial radius of the main halo. 

Weak gravitational lensing provides a powerful tool for 
studying the mass distribution of clusters of galaxies as well 
as the large-scale structure in the Universe (see Bartelmann 
& Schneider 2001; Refregier 2003; Schneider et al. 2006; 
Munshi et al. 2008, for reviews on weak lensing). Several 
methods were developed on cluster mass reconstructions 
(see e.g. Kaiser & Squires 1993; Seitz & Schneider 2001). 
Besides weak lensing, strong lensing systems are also em- 
ployed to extend the weak lensing analysis into the inner 
part of clusters, such as the multiply-imaged systems and 
the corresponding identification of the location of critical 
curves (Cacciato et al. 2006; Bradac et al. 2005b). Since 
then some cluster mass reconstructions have been success- 
fully carried out (Bradac et al. 2005a, 2006; Merten et al. 
2009). Both methods reconstruct the gravitational poten- 
tial if), based on minimizing a x 2 ~ function comparing ob- 
served data with model (Bartelmann 1996; Seitz et al. 
1998). 

Flexion has been recently studied as the derivative of the 
shear, and responds to small-scale variations in the gravi- 
tational potential (Goldberg & Natarajan 2002; Goldberg 
& Bacon 2005; Bacon et al. 2006). The two complex spin-1 
and spin-3 components of flexion together with the shear 



describe how intrinsically round sources are mapped into 
'arclets'. The measurement of flexion can be obtained in 
principle from the same imaging data as used for shear 
measurements, and methods for that have been discussed 
in the literature (Irwin & Shmakova 2006; Okura et al. 2007; 
Goldberg & Leonard 2007; Schneider & Er 2008). One dis- 
advantage of shear technique is that galaxies have a broad 
intrinsic ellipticity distribution which constitutes a serious 
noise component. Since third-order brightness moments of 
galaxies are suspected to be much smaller than second- 
order ones, we expect a lower intrinsic noise in flexion mea- 
surements. It has been noted that flexion can contribute 
to study the galaxy and cluster halo mass and density pro- 
files, and particularly sensitive to substructure (Bacon et al. 
2009; Leonard et al. 2009; Okura et al. 2008). Okura et al. 
(2008); Leonard et al. (2007) have used flexion measure- 
ments to detect substructure in Abell 1689. 

We propose a method that combines shear measure- 
ments and flexion into a weak lensing mass reconstruc- 
tion, which enable us to reconstruct the mass distribu- 
tion of a lens with special sensitivity to substructure. It is 
based on an estimations of the lensing potential using x 2 ~ 
minimization. The resulting equations are first linearzied 
and then solved by iteration. We outline the basic lens- 
ing formalism in Sect. 2, and the ideas of reconstruction 
method in Sect. 3 (calculations are shown in the Appendix). 
Numerical tests are presented in Sect. 4, and conclusions are 
given in Sect. 5. 

2. Basic formalism 

We mainly follow the notations of Schneider & Er (2008) 
to present the basic gravitational lensing shear and flexion 
formalism. We adopt the thin lens approximation where the 
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lensing mass distribution is projected onto the lens plane. 
Thus, all the lensing properties are contained in the lensing 
potential 0, which is the projected Newtonian gravitational 
potential on the lens plane 



i>(0) = - [ A 2 e'n{e')\n\e -e'\ , 



(1) 



where 6 denotes the (angular) position in the lens plane, 
and k(9) is the projected mass density E(0) of the lens in 
units of 



AirG D A D A . 



so 



K(0) = 



(2) 



for a fiducial source located at a redshift z — > oo. Here 
Dco,D A and -Dd.oo are the angular diameter distances be- 
tween the observer and the source, the observer and the 
lens and between the lens and the source, respectively. 
The second-order local expansion of the lensing equation in 
Cartesian coordinates reads /3j = 9i — — ip,ijkQjQk/2, 

where indices separated by a comma denote partial deriva- 
tives with respect to 9 i} and the Einstein summation con- 
vention is used. The surface mass density k and the com- 
plex shear 7 = 71 + i72 are given in terms of the deflection 
potential through 



k = 0,ii + V\22)/2, 7 = 0,11 - V>,22)/2 + 10. 



12 



(3) 



The two flexions J 7 , Q arc combinations of third-order 
derivatives of ip, and also related to the gradient of k and 

7 

1 



T = 

Q = 
T = 



[-0,111 + V\122 + i 0,112 + -0,222)] 

\ [V-,111 

v c k = v*7, g = v c 



3^,122 + i (3^,112 - 0,222)] 



=7, 



(4) 



where wc defined the differential operators 

„ d . d „„ d . d 
V r = h i ; V = i 

Mi d0 2 c d6i dd 2 

The second-order lensing equation in complex notation then 
reads 



P = (1 - k)0 - 7<9* - J T* 8 2 - \^T99* - (9*) 2 



(5) 



Due to the mass-sheet degeneracy in lensing system (see, 
e.g., Falco et al. 1985; Schneider & Seitz 1995) the lens 
equation is rewritten in reduced form as (see Schneider & 
Er 2008) 



(3 



and where \&i, ^3 are given by 
Gi - gGt 



4(1 

G3 
4 



99*) 
9(Gi 



gGt) 



Mi-gg*) 



(8) 



All the relation above are written for sources located at 
redshift z — > 00. For a source at redshift z and lens at 
redshift Zd, the lensing strength is reduced by multiplying ip 
(and all its derivatives) by a 'cosmological weight' function, 



Z(z) 



D, 



d.oo 



D_ 



H(z - z d ), 



(9) 



where D As and D s are the angular diameter distances be- 
tween the lens and the source, and the observer and the 
source. H(z — z A ) is the Heaviside step function to give 
zero weight to sources located in front of the lens. Thus, at 
redshift z 



k{z) 
T{z) 



Z{z)k, j(z) = 
Z(z)T, G(z) 



-~z(z)g. 



(io) 



Information about the lensing potential can be obtained 
through the reduced shear and reduced flexion by measur- 
ing the ellipticity e and the higher-order brightness mo- 
ments, or called HOLICs (Okura et al. 2007). The informa- 
tion on the lens potential is contained in the transformation 
between the source and image ellipticity 



e = < 



e-5 



9t* 



for \g\ < 1; 



for U > 1, 



(11) 



and the higher-order brightness moments 
9F -12Q 2 



Tj 3 w Ti 



T a - 



'-Gi, 



3F 



Go 



(12) 



where Q,T,F are the second-, third- and fourth- order 
brightness moments, respectively (see Schneider & Er 2008, 
for definitions), the superscript s indicates source quanti- 
ties, and the subscript indicates the spin of the moments. 
If we assume that the ensemble average of intrinsic mo- 
ments vanish, (e s ) = 0, (Tf 3 ) = 0, the expectation value 
for the galaxy image quantities yields the estimator of re- 
duced shear and reduced flexion 



(l-«) 
9 - qO* 



*I(Gi) 9 l - 2<F 1 (G 1 ) 99* - * 3 (Gi,G 3 ) (9*) 



(6) 
2 



where g, G\,Ga are the reduced shear and the reduced flex- 
ion components, 
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G l = 
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G 1 



G, 



9F - 12Q 2 
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» if|s|<l, 
(l/e*), otherwise, 



(13) 



3F " 



-Ta 



(14) 



F + gF* 
(l-«) 

g + gj- 



(7) 



Here we use ti.a to indicate the estimators of the reduced 
flexion. These estimates of reduced shear and reduced flex- 
ion break down when the image is located on or close to a 
critical curve. 
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3. Mass reconstruction method 

The strategy of our method is similar to that of the strong- 
and weak-lensing united method (Bradac et al. 2005b, a). 
Thus, we include flexion in the reconstruction of the lens- 
ing potential ip. There are two reasons for combining shear 
and flexion other than using flexion alone in the mass re- 
construction. In current observations, the number density 
of images from which flexion can be estimated is low, as 
we can see in our simulation data. Even with future higher 
flexion number density surveys, flexion is most sensitive 
to substructures, but rather insensitive to the larger-scale 
(cluster-scale) smooth potential. 

3.1. The x 2 -Unction 

We describe the cluster mass distribution by the deflection 
potential ip on a regular grid, use finite differencing to calcu- 
late the deflection angle, the reduced shear and the flexion. 
These quantities on the grid are considered to describe the 
lens model which can be compared with data. Our aim is 
to seek a potential ip that minimizes the difference between 
model (g, Gi, G 3 ) and data (e, t\,ts). We therefore define a 
% 2 -function 

x 2 (V0 = x 2 (V0 + XeW) + x/WO + vR{>l>), (is) 

where x 2 : Xe(V0 an d X/(V0 contain the information from 
strong lensing multiple image systems, shear and flexion, 
respectively. r/R is a regularization term, used to smooth 
out small-scale noise components, i.e., to avoid overfitting 
the data. For a multiple image system with N m images 
located at 6 m , the x 2 term is well defined by Bradac et al. 
(2005b) 

N m 

.2 _ V- V- uT a-lu (16) 



X, S 



J2 E ^b mi 



str m— 1 



where X^str i s the sum over all strong lens system, and 
b m = S m — a(0 m ) — (3 S . (3 S is the average source position, 
and S is the covariance matrix S = diag(cr 2 1 , cr 2 2 )- The 
shear term \\ is 



E 

i=i 



\e t -g(e, Zl )f 



(17) 



where N g is the number of background galaxies, and 



a 2 = (l - \g\ 2 ) 2 a% + a^, (18) 



with a> w 0.3 is the standard deviation (Brainerd et al. 
1996), er orr is the measurement error, which we take 0.1 in 
this paper, and g is the reconstructed reduced shear. 

The flexion term is defined in a way similar to the shear 



AT 

g/ l^-G^,^)! 2 | |%-G 3 ^,z,)| 2 j 
i=i ^ 



<4 



where Gi,3 are the reduced flexions at position Oi and red- 
shift Zi in our model. Note that we can measure only the 
reduced flexion. There is a significant difference between the 
reduced flexion and its approximation J 7 / (1 — k) or G/{1— n) 
in the region where the shear is not small. However, the 
definition of reduced flexion (Eq.7) renders the x 2 function 



(Eq.19) complicated and the equations become difficult to 
solve. We thus define 



1-k' 



G' 



1 - K 



From Eq.(7) it is easy to obtain 



G i = TT^ Gi=G 3 -gG[. 



(20) 



(21) 



We thus use a modified estimator for the observed reduced 
flexion, 



, _ h - et\ _ 



(22) 



and replace t\, t 3 by t[, t' 3 in (Eq.19). The flexion term \ 2 
is thus redefined as 



N, 



( 



i=l 



+ 



'tv 



't3> 



V 



(23) 



The al v and cr 2 3 / are different from the dispersion in flexion 
measurement or intrinsic flexion variance, which are both 
difficult to obtain from current observations. In Goldberg 
Leonard (2007), the intrinsic scatter of flexion are esti- 



mated as o a \j r \ — 0.03 and a, 



a\Q\ 



0.04, where a is the 



semi-major axis of the lensed image, so that the combina- 
tion a\T\ represents a dimensionless term. Hawken & Bridle 
(2009); Leonard et al. (2009) used a conservative estimate 
(jjr = 0.1/". Here we use 



'tv 



't3' 



dG[ 



dg 
dG'* 



dg 



a e + <T t3' 



(24) 



where an, Otz are the dispersion of flexion data. The first 
term on the right is scatter due to intrinsic ellipticity and 
shear measurement error. 

In principle, the total number of flexion measurements 
is the same as the number of galaxy images Nf = N g from 
which shear can be estimated. But in reality, flexion is mea- 
surable only for a subset of (larger) images; beside that, 
we discard low signal-to-noise flexion measurements. Thus 
the number of images from which flexion can be estimated 
is usually smaller than the number of galaxy images with 
shear measurements. 

To find the minimum x 2 -function, we solve the set of 
equations 

(25) 

which is in general a non-linear set of equations. This 
problem is solved by an iterative procedure (Bradac et al. 
2005b). We perform a three- level iteration process, start 
with initial model the fixed regularization, linearize 

the system (Appendix) , solve the linear system of equations 
as the inner-level. These steps are repeated until conver- 
gence of k is achieved. The middle-level is repeating the 
inner level with new regularization, based on the estimate 
of k from previous step result, until we can get Xrcd ~ 1- 
Finally we increase the number of grid points in the field 
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and repeat the first two level iterations on the new grid until 
we reach the final desired grid spacing, i.e., resolution. 
For the regularization we choose 



R 



N X ,Ny 

£ 



(n) 

ij 



(26) 



where k'") and k^" -1 ) are current and previous resulting k 
map in every middle- level. For the first step, we can use an 
initial model k*- ' which is obtained from other methods, or 
simple set = across the whole field. 

3,2. Grid point field 

Finite differencing techniques (Abramowitz & Stegun f 972) 
provide a way to calculate k and 7 on the potential ip grid 
field. We need 9 grid points for second-order derivatives 
required for k and 7, and 16 points for flexion. But by using 
4x4 points, one would obtain the flexion value at the center 
of a grid cell. Thus we instead use 5x5 finite difference 
scheme to obtain the flexion on the grid. For instance, 



1 



+2 

A it A,=-2 



WA < ,A J ^(* + Ai,j + A i ), (27) 



where the coefficients Wa 4 ,Aj used for spin-1 flexion are 
given in the left panel of Fig. 1. Bilinear interpolation can 
be used to obtain the shear and flexion for an arbitrary 
position in the field. We need to extend two grid rows and 
columns on each boundary of the entire field, which means 
we use a (TV + 4) 2 grid to perform the method for inner N 2 
grid. Note that the flexion drops quickly with increasing 
distance from the cluster center, so that there is few images 
with flexion signal near the grid boundary. 

Since shear and flexion are second- and third-order 
derivatives of ip, the potential field is not fixed with shear 
and flexion constraints alone. However, k is unchanged 
under the transformation ip(9) — > ip(0) + ipo + a. ■ 0, 
where ipo an d a are arbitrary constants. We leave the con- 
stants free to simplify the numerical solution. The mass- 
sheet degeneracy transformation of the potential is given 
by ip — > ip' = (1 — X)9 2 /2 + Xip, which does not affect the 
reduced shear and reduced flexion, but affects k as 



(1-A)+Ak, 



(28) 



which we call the A transformation. We will take advantage 
of this transformation for adjusting k later. We note that 
the mass-sheet transformation is exact only for all sources 
being located at the same redshift; for a redshift distribu- 
tion, the mass-sheet degeneracy is broken. However, this 
breaking is not very strong, as long as no multiple images 
of sources with very different redshifts are considered (see 
Bradac et al. 2005b) 

4. Numerical test of field properties 

4.1. NIS model 

First we perform our method on a Non-singular Isothermal 
Sphere model. The data is generated on a 40 x 40 grid 
potential field with 



tf>{6) 



(29) 



n 



($) — (z$)~ ( 



LI 



J_J 




T\ ->■ 1/(12A 3 ) 

Fig. 1. The finite differencing coefficients of .Fi(left) and 
CJi(right). The coefficients of Ti^Q-i) are the same as that of 
3~\(Q\) after rotating by 7r/2 anticlockwise(clockwise). 



where 0-^ = 1 arcmin and C = 0.2 arcmin. The reduced 
shear and reduced flexion are calculated by finite differenc- 
ing. To ensure that no strongly lensed objects are included, 
we reject all the data points of which the absolute value 
of the reduced shear is larger than 0.9. Fig. 2 shows the 
original and reconstructed radial k and \g\ profiles. The 
value is obtained by annular bin average on the grid field, 
the origianl one as well as the reconstructed one. The small 
fluctuations of the results around the input profile is mainly 
due to noise and low grid resolution. In the left panel, we 
plot the convergence k after applying the A transformation 
Eq.28. We can see that convergence is well recovered for 
large 9, up to the vertical line. For smaller 6, the results 
from shear alone deviate strongly from the input. However, 
using shear and flexion together significantly improves the 
resulting mass profile. In the right panel, the absolute value 
of the reduced shear is shown, which is not affected by mass- 
sheet degeneracy. Both results, from shear alone and from 
combining shear and flexion show agreement with the in- 
put model. Again for small 8, the result from the combined 
shear and flexion data gives better agreement with the in- 
put line. 

4.2. Mock data 

Our mock data uses a cluster taken from N-body simula- 
tions by Jing & Suto (2002); Jing (2002). The cluster is sim- 
ulated in the framework of a ACDM model with = 0.7, 
fl nl = 0.3, the normalization of power spectrum <7s = 0.9 
and Hubble constant Hq = 70kms _1 Mpc _1 . Dark matter 
halos are identified with the friends-of-friends method us- 
ing a linking length equal to 0.2 times the mean particle 
separation. The halo mass M is defined as the mass en- 
closed within the virial radius according to the spherical 
collapse model (Kitayama & Suto 1996; Bryan & Norman 
1998; Jing & Suto 2002). The virial mass of the cluster 
which we use here is 3.4 x 10 14 1i _1 Mq and its redshift is 
Zd = 0.326. The particles within a box with side length of 
2 virial radii were projected onto the lens plan. The sur- 
face densities are calculated using the smoothed particle 
hydrodynamics smoothing algorithm (Monaghan 1992) on 
a 4096 x 4096 grid. The lensing potential is obtained us- 
ing the fast Fourier transform method (Bartelmann et al. 
1998). 

A finer grid 0.1 arcsec resolution is created around the 
cluster center with side length of 4 arcmin. We obtain the 
mapping (3(0) from the lens equation for each grid point 
and the corresponding second-order derivatives using cubic 
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0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 

e(arcmin) 6(arcmin) 

Fig. 2. Radial profile of the NIS cluster, the solid line is the input model, the plus points are the reconstructed result with weak 
lensing shear only, and the cross points are the result with weak lensing shear and flexion. The mock data is generated out of the 
vertical line. Left panel: convergence k after A transformation, Right panel: the absolute value of reduced shear \g\, all for a source 
at z s — > oo. 



spline interpolation. We perform cubic spline interpolation 
again to the second-order derivatives to get the third-order 
derivatives. The background sources are located at random 
position in the source planes, and their redshifts follow a 
Gamma distribution 

z 2 

p( z ) = ex P (- z / z o) , (30) 
z z 

where zq = 1/3. The peak is at z = 2/3 and the mean 
rcdshift is (z) = 3zo = 1- The Newton- Raphson method 
is used to find the corresponding images position on the 
image plane. The reduced shear and flexion for each image 
were linear interpolated using four nearest grid points. 

One simulated cluster with two different projection di- 
rections is used to generate the two different sets of mock 
data (Fig. 3). We name the two sets of data as d02 and 
d03. All of the reduced shear data is used in our calcula- 
tion, but for flexion, we only consider reduced flexion with 
absolute value in the range from 0.01 to 0.5. For higher 
values, merging of multiple images most likely render any 
flexion measurement impossible (in fact, there the whole 
concept of flexion break down - see Schneider & Er 2008). 
On the other hand, very small values of the reduced flexion 
are exceedingly difficult to measure and would contribute 
very little information due to the low signal-to-noise. 

4.3. Reconstructed k map 

The two mock catalogues are used to test the performance 
of our method. We start with an initial 20 x 20 grid, in- 
crease N x and N y by one for each outer loop iteration, up 
to a 30 x 30 grid. We use N g = 1000 weak lensing galaxy 
images in each reconstruction, which is an accessible back- 
ground galaxy number density of ~ 60 images arcmin~ 2 . 
The result of the reconstructions are shown in Fig. 4 for 
d02 and Fig. 5 for d03. The initial regularization parameter 
is set to i] = 200 for cluster d02, and r\ — 300 for cluster 
d03, and is increased by a factor of 10 for each outer-level 



iteration. It is usually better to set high rj and allowing k to 
change slowly. Since our reconstruction is done in a three- 
level iteration, and in each step we ensure x 2 /-^dof ~ 1j the 
method can successfully adapt to the data and the results 
arc not sensitive to the initial choice of r\. We need also 
an initial k° field for the regularization. A simple model 
k° = 0.01 is used here. We have also performed additional 
reconstructions with different initial models, and found the 
results are nearly independent of the initial k°, but a more 
realistic model allows for a faster convergence. 

The number of flexion values that we used for the two 
cluster is different. We obtained 89 reduced flexion value 
between 0.01 to 0.5 in d02 (Fig. 4), and 60 that for d03 
(Fig. 5). Since there is more significant extended structure 
in cluster d02 than in d03, we found that there are more 
high flexion signal data (the absolute value of the spin-1 
reduced flexion greater than 0.01) in cluster d02 than that 
in cluster d03. 

The results show that our method can reproduce the 
main properties of the projected mass distribution of both 
clusters, and is especially powerful in resolving the sub- 
structure. Fig. 4 shows that our method can reconstruct the 
k map by combining multiple images, weak lensing shear 
and flexion data for cluster d02. Unfortunately we cannot 
clearly distinguish local small-scale maxima which are due 
to noise from the true low-mass substructure, even with the 
help of flexion, like the one in the bottom corner of Fig. 4. 
However we can sec that after including flexion, the halos 
become more peaky and more massive substructures are 
resolved and at the correct positions. In Fig. 5 it is encour- 
aging to see that besides the overall mass distribution of 
the cluster halo, our method can resolve the small clumps 
after including flexion data. The small clump is not very 
significant even in the input convergence (Fig. 3), which is 
at a resolution of 400 x 400. We also performed additional 
tests in which we use different sets of weak lensing shear 
and flexion data for both clusters and confirmed the valid- 
ity of our method. In some cases of cluster d03 data, the 
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Fig. 3. The convergence map of two simulated cluster used for generating mock strong lensing multiple images, weak lensing shear 
and flexion data. The triangles denote the image positions of multiple image systems which we use for the reconstruction, n is 
plotted for sources at z s — > oo, given in linear gray-scale and contours. We name the left panel cluster d02 and the right one cluster 
d03. 
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Fig. 4. K-maps obtained from shear and flexion reconstruction of the mock data cluster d02 after the A transformation (Eq.28), 
where A is chosen to yield the smallest D 2 (Eq.31). The left panel shows the result using JV g = 1000 shear galaxies and 3 multiple 
image systems, while in the right panel we add 89 flexion values. 



shape properties can be better reconstructed and the small 
clump can be clearly resolve, but the position of the small 
clump can have an offset up to 10 arcsec from its real po- 
sition. In some other realizations, the small clump cannot 
be clearly resolved, which is due to the noise and local low 
background image density around the small clump. 

As an additional test, we calculate the difference D 2 be- 
tween our reconstruction and the input k, which is defined 

as 

^ 2 = ^f>,-4 nput) i 2 > (31) 



where N' is number of grid points which are not inside 
the critical curves. We apply the A transformation (Eq.28) 
to the convergence results, choosing A such as to yield the 
smallest difference D 2 . For cluster d02, the D 2 W of the re- 
sult from shear only is 0.0086, and when including flex- 
ion, it becomes 0.0085. For cluster d03, D 2 W = 0.0077 and 
Dg wf = 0.0066, respectively. Whereas we have employed 
the multiple image systems in our reconstruction, in cases 
where such strong lensing systems are not available the rel- 
ative improvement obtained by the inclusion of flexion is 
expected to be considerable larger. 
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Fig. 5. Reconstructed K-maps from mock data cluster d03 after the A transformation. Left panel shows the result using iVg = 1000 
shear galaxies, and 3 multiple image systems, whereas the right panel also uses 60 flexion values. 



We also enlarged the allowed range for flexion to 
[0.001,0.5] and find that in some cases, D 2 becomes larger 
after including flexion. After checking our resulting k on 
all grids, we identify the points which give larger D 2 after 
combining the flexion signal. We are aware that the rela- 
tion of the reduced flexion and the reduced shear (Eq.22) 
can transfer extra noise from shear into flexion, especially in 
case of large intrinsic shear noise. In brief, if we decrease the 
threshold of flexion and use more data with lower signal-to- 
noise. The minimizing of x 2 will relatively lose the constrain 
on the region with strong flexion signal, thus a bigger D 2 
is obtained. 

Finally a word on the dispersion of flexion a' t . This is a 
difficult parameter to determine at the moment, since we 
have little knowledge about the noise behavior in flexion 
measurement. The one we used in this paper (Eq.24) has a 
problem: as pointed out by Bacon et al. (2009), the flexion 
variance is biased by the content of substructure. Moreover, 
we noticed that the oy we used is underestimated, it can 
be seen from following: \ 2 /Nf is significant smaller than 
1, i.e. at' is too large and the constrain from flexion is lose 
and down weighted. If we apply more steps of iteration, the 
noise might be over fitted. In that case, the cluster becomes 
very peaky, and looks like being truncated at some edge 
region. Some other forms for a' t have been also tested, e.g. 
in analogy to the shear 



al = {l-{6 \t'\f)\l+a: 



(32) 



where 9q is size of the image. We can easily see that this 
Uf is not independent of image size and a' crr is different 
from that of shear (Eq.18), since it is not dimcnsionless. 
It turns out that flexion is underestimated, of which the 
parameter we used is a t s = 0.1/" and a' m = 0.1/". Thus, 
flexion variance will affect the constraint of flexion. A more 
suitable model of flexion noise need to be constructed. 



5. Conclusions 

In this paper we propose a method for projected cluster 
mass reconstruction, which combines strong, weak lens- 
ing shear and flexion data. The method is based on a 
least-x 2 fitting of the lensing potential ip (Bradac et al. 
2005b; Cacciato et al. 2006). The particular strength of this 
method is that the flexion data provides more information 
to the inner parts of the cluster and on substructure. 

We test the performance of the method on our mock 
clusters, comparing the results with and without flexion. In 
the NIS cluster, our method can reproduce the radial pro- 
file 1 of the convergence and the reduced shear. Flexion can 
significantly improve the result in the inner part of the mass 
profile. In the other test, we generate our mock data from 
simulated clusters. We are able to reconstruct the main 
properties of the cluster mass distributions; in particular 
when the flexion data is included, our method can success- 
fully resolve the cluster and substructure. In addition, our 
result is almost independent of the initial model k° and the 
regularization parameter r\. 

We have assumed that the intrinsic flexion is small. 
However, the correction for the reduced flexion introduces 
extra noise from shear to flexion, especially in the case 
where the intrinsic galaxies are highly elliptical. The ef- 
fect of the noise is limited where the shear and flexion are 
strong, which is, however, not the case in the outer regions 
of the clusters. It is of interest to study the relation of in- 
trinsic noise and flexion in detail. 

In Leonard et al. (2007); Okura et al. (2008), the result 
of mass reconstruction by flexion has shown that flexion 
is sensitive to substructure, and insensitive to the smooth 
component of the cluster. Our method of combining shear 
and flexion takes the advantages of shear on cluster scale, 
and of flexion for substructures. 

The number density of background galaxies that we 
used, <~ 60 images arcmin" 2 , is accessible by current tele- 
scopes such as HST. With future telescopes, the accuracy 
and number density of images can be improved. Once flex- 
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ion can be measured accurately, the noise behavior and in- 
trinsic flexion scatter are certainly need to be studied before 
putting flexion into practical use. 
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Appendix A: The iteration 

We present the method outlined in Sect. 3.1 on how we lin- 
earize and solve Eq.(25). The lcnsing quantities are calcu- 
lated by finite differencing, and are thus linear combinations 
of ip at each position. They are expressed in the following 
matrix notations (Fig. 1) 



7i(^) = M7^fe; 



K{6i)=M^ k ; 
72^) = M^Vfe! 



(A.l) 
(A.2) 

(A.3) 

(A.4) 



Then we plug these into Eq.(25) and obtain the full form 
of the equations. Here we write down the first flexion term 
as an example. The strong lens multiple images, shear and 
regularization part can be found in Bradac et al. (2005b) 
and the second flexion term will be analogously to the first 
one, 

Xfi WO = E M_.^2-2 i ( A - 5 ) 
»=1 ^ 



'tv 



here t\ is given by Eq.(22), and we drop the prime on t\ in 
this appendix, for notational simplicity. We also omit the 
index i to all parameters of every galaxy for simplicity. We 
fix the denominator a\ x = (1 — «0 2 ofi as constant at each 
iteration step, so they will not appear in the derivative, 



dipk 



Ni 

i=i 

N f 



tl r — 



E 

1,2 



[(1 - K)t lr - F r ] h 



dn 

dipt 



with the matrix Bkj and vector Vj containing the contri- 
butions from the nonlinear part. 



K = — [ M t] M ll 



Ml] Mti 



tniMnMg + MftMZ) 



+ t 12 {M%M? k + M%M$) + (i?! + t? 2 )M£.M&] , (A.8) 

where i denote summation over all the galaxy images. The 
data vector is 



Vf = E i [*n^y + h2M% + (t 2 u + t\ 2 )M^ . (A.9) 
<=i CT « 

Here B k - and Vf are contributions from spin-1 flexion. 
The same calculation can be performed for spin-3 flexion, 
and a similar result is obtained. Combining all the matrixes 
and data vectors, one can complete the matrix Bkj and the 
vector Vj. 



Ok dF r ^ 



i=i 



M]\M, 



+ Mlj Ml 



^^ + (*?i+*?2)M5M^™-Tval. 
+t 11 {Mi]M? k + MtlM? j ) 



+t 12 (M t i jM? k 



a 



ii 
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